
<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta charset="utf-8" />
    <title>lmp &#8212; Texas A&amp;M Oil spill / Outfall Calculator 2.1.0 documentation</title>
    <link rel="stylesheet" href="../_static/alabaster.css" type="text/css" />
    <link rel="stylesheet" href="../_static/pygments.css" type="text/css" />
    <script id="documentation_options" data-url_root="../" src="../_static/documentation_options.js"></script>
    <script src="../_static/jquery.js"></script>
    <script src="../_static/underscore.js"></script>
    <script src="../_static/doctools.js"></script>
    <script src="../_static/language_data.js"></script>
    <link rel="index" title="Index" href="../genindex.html" />
    <link rel="search" title="Search" href="../search.html" />
   
  <link rel="stylesheet" href="../_static/custom.css" type="text/css" />
  
  
  <meta name="viewport" content="width=device-width, initial-scale=0.9, maximum-scale=0.9" />

  </head><body>
  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          

          <div class="body" role="main">
            
  <h1>Source code for lmp</h1><div class="highlight"><pre>
<span></span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Lagrangian Multiphase Plume</span>
<span class="sd">===========================</span>

<span class="sd">This module contains the numerical solution for the `bent_plume_model` module.</span>
<span class="sd">Some of the general tools for handling the multiphase components are</span>
<span class="sd">contained in `dispersed_phases`.</span>

<span class="sd">&quot;&quot;&quot;</span>
<span class="c1"># S. Socolofsky, November 2014, Texas A&amp;M University &lt;socolofs@tamu.edu&gt;.</span>
<span class="kn">from</span> <span class="nn">__future__</span> <span class="kn">import</span> <span class="p">(</span><span class="n">absolute_import</span><span class="p">,</span> <span class="n">division</span><span class="p">,</span> <span class="n">print_function</span><span class="p">)</span>

<span class="kn">from</span> <span class="nn">tamoc</span> <span class="kn">import</span> <span class="n">seawater</span>
<span class="kn">from</span> <span class="nn">tamoc</span> <span class="kn">import</span> <span class="n">dispersed_phases</span>

<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">scipy</span> <span class="kn">import</span> <span class="n">integrate</span>
<span class="kn">from</span> <span class="nn">copy</span> <span class="kn">import</span> <span class="n">deepcopy</span>

<div class="viewcode-block" id="derivs"><a class="viewcode-back" href="../autodoc/bpm/lmp.derivs.html#lmp.derivs">[docs]</a><span class="k">def</span> <span class="nf">derivs</span><span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">q</span><span class="p">,</span> <span class="n">q0_local</span><span class="p">,</span> <span class="n">q1_local</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">particles</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Calculate the derivatives for the system of ODEs for a Lagrangian plume</span>
<span class="sd">    </span>
<span class="sd">    Calculates the right-hand-side of the system of ODEs for a Lagrangian </span>
<span class="sd">    plume integral model.  The continuous phase model matches very closely </span>
<span class="sd">    the model of Lee and Cheung (1990), with adaptations for the shear </span>
<span class="sd">    entrainment following Jirka (2004).  Multiphase extensions following the</span>
<span class="sd">    strategy in Socolofsky et al. (2008) with adaptation to Lagrangian plume</span>
<span class="sd">    models by Johansen (2000, 2003) and Yapa and Zheng (1997).  </span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    t : float</span>
<span class="sd">        Current value for the independent variable (time in s).</span>
<span class="sd">    q : ndarray</span>
<span class="sd">        Current value for the plume state space vector.</span>
<span class="sd">    q0_local : `bent_plume_model.LagElement`</span>
<span class="sd">        Object containing the numerical solution at the previous time step</span>
<span class="sd">    q1_local : `bent_plume_model.LagElement`</span>
<span class="sd">        Object containing the numerical solution at the current time step</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the simulation.</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the bent</span>
<span class="sd">        plume model.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `bent_plume_model.Particle` objects containing the dispersed </span>
<span class="sd">        phase local conditions and behavior.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    yp : ndarray</span>
<span class="sd">        A vector of the derivatives of the plume state space.</span>
<span class="sd">    </span>
<span class="sd">    See Also</span>
<span class="sd">    --------</span>
<span class="sd">    calculate</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Set up the output from the function to have the right size and shape</span>
    <span class="n">qp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">q</span><span class="o">.</span><span class="n">shape</span><span class="p">)</span>
    
    <span class="c1"># Update the local Lagrangian element properties</span>
    <span class="n">q1_local</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">q</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">particles</span><span class="p">)</span>
    
    <span class="c1"># Get the entrainment flux</span>
    <span class="n">md</span> <span class="o">=</span> <span class="n">entrainment</span><span class="p">(</span><span class="n">q0_local</span><span class="p">,</span> <span class="n">q1_local</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
    
    <span class="c1"># Get the dispersed phase tracking variables</span>
    <span class="p">(</span><span class="n">fe</span><span class="p">,</span> <span class="n">up</span><span class="p">,</span> <span class="n">dtp_dt</span><span class="p">)</span> <span class="o">=</span> <span class="n">track_particles</span><span class="p">(</span><span class="n">q0_local</span><span class="p">,</span> <span class="n">q1_local</span><span class="p">,</span> <span class="n">md</span><span class="p">,</span> <span class="n">particles</span><span class="p">)</span>
    
    <span class="c1"># Conservation of Mass</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">md</span>
    
    <span class="c1"># Conservation of salt and heat</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">md</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">Sa</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">md</span> <span class="o">*</span> <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">()</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">Ta</span>
    
    <span class="c1"># Conservation of continuous phase momentum.  Note that z is positive</span>
    <span class="c1"># down (depth).</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">md</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">ua</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="n">md</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">va</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span> <span class="n">p</span><span class="o">.</span><span class="n">g</span> <span class="o">/</span> <span class="p">(</span><span class="n">p</span><span class="o">.</span><span class="n">gamma</span> <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">rho_r</span><span class="p">)</span> <span class="o">*</span> <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">Fb</span> <span class="o">+</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">M</span> <span class="o">*</span> 
            <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">-</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">rho</span><span class="p">))</span> <span class="o">+</span> <span class="n">md</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">wa</span>
    
    <span class="c1"># Constant h/V thickeness to velocity ratio</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">6</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
    
    <span class="c1"># Lagrangian plume element advection (x, y, z) and s along the centerline</span>
    <span class="c1"># trajectory</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">7</span><span class="p">]</span> <span class="o">=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">u</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">8</span><span class="p">]</span> <span class="o">=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">v</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">9</span><span class="p">]</span> <span class="o">=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">w</span>
    <span class="n">qp</span><span class="p">[</span><span class="mi">10</span><span class="p">]</span> <span class="o">=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">V</span>
    
    <span class="c1"># Conservation equations for each dispersed phase</span>
    <span class="n">idx</span> <span class="o">=</span> <span class="mi">11</span>
    
    <span class="c1"># Track the mass dissolving into the continuous phase per unit time</span>
    <span class="n">dm</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">nchems</span><span class="p">)</span>
    
    <span class="c1"># Compute mass and heat transfer for each particle</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">)):</span>

        <span class="c1"># Only simulate particles inside the plume</span>
        <span class="k">if</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">integrate</span><span class="p">:</span>
            
            <span class="c1"># Dissolution and Biodegradation</span>
            <span class="k">if</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">issoluble</span><span class="p">:</span>
                <span class="c1"># Dissolution mass transfer for each particle component</span>
                <span class="n">dm_pc</span> <span class="o">=</span> <span class="o">-</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">A</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nbe</span> <span class="o">*</span> \
                           <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">beta</span> <span class="o">*</span> <span class="p">(</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">Cs</span> <span class="o">-</span> 
                           <span class="n">q1_local</span><span class="o">.</span><span class="n">c_chems</span><span class="p">)</span> <span class="o">*</span> <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
                
                <span class="c1"># Update continuous phase temperature with heat of </span>
                <span class="c1"># solution</span>
                <span class="n">qp</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">+=</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">dm_pc</span> <span class="o">*</span> \
                         <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">neg_dH_solR</span> \
                         <span class="o">*</span> <span class="n">p</span><span class="o">.</span><span class="n">Ru</span> <span class="o">/</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">M</span><span class="p">)</span>
                
                <span class="c1"># Biodegradation for for each particle component</span>
                <span class="n">dm_pb</span> <span class="o">=</span> <span class="o">-</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">k_bio</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">m</span> <span class="o">*</span> \
                     <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nbe</span> <span class="o">*</span> <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
                
                <span class="c1"># Conservation of mass for dissolution and biodegradation</span>
                <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">:</span><span class="n">idx</span><span class="o">+</span><span class="n">q1_local</span><span class="o">.</span><span class="n">nchems</span><span class="p">]</span> <span class="o">=</span> <span class="n">dm_pc</span> <span class="o">+</span> <span class="n">dm_pb</span>
                
                <span class="c1"># Update position in state space </span>
                <span class="n">idx</span> <span class="o">+=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">nchems</span>
                
            <span class="k">else</span><span class="p">:</span>
                <span class="c1"># No dissolution</span>
                <span class="n">dm_pc</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">nchems</span><span class="p">)</span>
                <span class="c1"># Biodegradation for insoluble particles</span>
                <span class="n">dm_pb</span> <span class="o">=</span> <span class="o">-</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">k_bio</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">m</span> <span class="o">*</span> \
                    <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nbe</span> <span class="o">*</span> <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
                <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="n">dm_pb</span>
                <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
            
            <span class="c1"># Update the total mass dissolved</span>
            <span class="n">dm</span> <span class="o">+=</span> <span class="n">dm_pc</span>
            
            <span class="c1"># Heat transfer between the particle and the ambient</span>
            <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">A</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nbe</span> <span class="o">*</span> \
                        <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">rho_p</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">cp</span> <span class="o">*</span> \
                        <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">beta_T</span> <span class="o">*</span> <span class="p">(</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">T</span> <span class="o">-</span> \
                        <span class="n">q1_local</span><span class="o">.</span><span class="n">T</span><span class="p">)</span> <span class="o">*</span> <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
            
            <span class="c1"># Heat loss due to mass loss</span>
            <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">+=</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">dm_pc</span> <span class="o">+</span> <span class="n">dm_pb</span><span class="p">)</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">cp</span> <span class="o">*</span> \
                       <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">T</span>
            
            <span class="c1"># Take the heat leaving the particle and put it in the continuous </span>
            <span class="c1"># phase fluid</span>
            <span class="n">qp</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">-=</span> <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
            
            <span class="c1"># Particle age</span>
            <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
            
            <span class="c1"># Follow the particles in the local coordinate system (l,n,m) </span>
            <span class="c1"># relative to the plume centerline</span>
            <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
            <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">up</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">fe</span> <span class="o">*</span> <span class="n">q</span><span class="p">[</span><span class="n">idx</span><span class="p">])</span> <span class="o">*</span> <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
            <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">up</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">-</span> <span class="n">fe</span> <span class="o">*</span> <span class="n">q</span><span class="p">[</span><span class="n">idx</span><span class="p">])</span> <span class="o">*</span> <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span>
            
        <span class="k">else</span><span class="p">:</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">nc</span> <span class="o">+</span> <span class="mi">5</span>
    
    <span class="c1"># Conservation equations for the dissolved constituents in the plume</span>
    <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">:</span><span class="n">idx</span><span class="o">+</span><span class="n">q1_local</span><span class="o">.</span><span class="n">nchems</span><span class="p">]</span> <span class="o">=</span> <span class="n">md</span> <span class="o">/</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">ca_chems</span> \
        <span class="o">-</span> <span class="n">dm</span> <span class="o">-</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">k_bio</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">cpe</span>
    <span class="n">idx</span> <span class="o">+=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">nchems</span>
    
    <span class="c1"># Conservation equation for the passive tracers in the plume</span>
    <span class="n">qp</span><span class="p">[</span><span class="n">idx</span><span class="p">:]</span> <span class="o">=</span> <span class="n">md</span> <span class="o">/</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">ca_tracers</span>
    
    <span class="c1"># Return the slopes</span>
    <span class="k">return</span> <span class="n">qp</span></div>


<div class="viewcode-block" id="calculate"><a class="viewcode-back" href="../autodoc/bpm/lmp.calculate.html#lmp.calculate">[docs]</a><span class="k">def</span> <span class="nf">calculate</span><span class="p">(</span><span class="n">t0</span><span class="p">,</span> <span class="n">q0</span><span class="p">,</span> <span class="n">q0_local</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">derivs</span><span class="p">,</span> <span class="n">dt_max</span><span class="p">,</span> 
              <span class="n">sd_max</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Integrate an the Lagrangian plume solution</span>
<span class="sd">    </span>
<span class="sd">    Compute the solution tracking along the centerline of the plume until </span>
<span class="sd">    the plume reaches the water surface, reaches a neutral buoyancy level </span>
<span class="sd">    within the intrusion layer, or propagates a given maximum number of</span>
<span class="sd">    nozzle diameters downstream.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    t0 : float</span>
<span class="sd">        Initial time (s)</span>
<span class="sd">    q0 : ndarray</span>
<span class="sd">        Initial values of the state space vector</span>
<span class="sd">    q0_local : `bent_plume_model.LagElement`</span>
<span class="sd">        Object containing the numerical solution at the initial condition</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the simulation.</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the bent</span>
<span class="sd">        plume model.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `bent_plume_model.Particle` objects containing the dispersed </span>
<span class="sd">        phase local conditions and behavior.</span>
<span class="sd">    derivs : function handle</span>
<span class="sd">        Pointer to the function where the derivatives of the ODE system are</span>
<span class="sd">        stored.  Should be `lmp.derivs`.</span>
<span class="sd">    dt_max : float</span>
<span class="sd">        Maximum step size to use in the simulation (s).  The ODE solver </span>
<span class="sd">        in `calculate` is set up with adaptive step size integration, so </span>
<span class="sd">        this value determines the largest step size in the output data, but </span>
<span class="sd">        not the numerical stability of the calculation.</span>
<span class="sd">    sd_max : float</span>
<span class="sd">        Maximum number of nozzle diameters to compute along the plume </span>
<span class="sd">        centerline (s/D)_max.  This is the only stop criteria that is user-</span>
<span class="sd">        selectable.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    t : ndarray</span>
<span class="sd">        Vector of times when the plume solution is obtained (s).</span>
<span class="sd">    y : ndarray</span>
<span class="sd">        Matrix of the plume state space solutions.  Each row corresponds to</span>
<span class="sd">        a time in `t`.</span>
<span class="sd">    </span>
<span class="sd">    See Also</span>
<span class="sd">    --------</span>
<span class="sd">    derivs, bent_plume_mode.Model</span>
<span class="sd">        </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Create an integrator object:  use &quot;vode&quot; with &quot;backward </span>
    <span class="c1"># differentitaion formula&quot; for stiff ODEs</span>
    <span class="n">r</span> <span class="o">=</span> <span class="n">integrate</span><span class="o">.</span><span class="n">ode</span><span class="p">(</span><span class="n">derivs</span><span class="p">)</span><span class="o">.</span><span class="n">set_integrator</span><span class="p">(</span><span class="s1">&#39;vode&#39;</span><span class="p">,</span> <span class="n">method</span><span class="o">=</span><span class="s1">&#39;bdf&#39;</span><span class="p">,</span> <span class="n">atol</span><span class="o">=</span><span class="mf">1.e-6</span><span class="p">,</span>
        <span class="n">rtol</span><span class="o">=</span><span class="mf">1e-3</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="mi">5</span><span class="p">,</span> <span class="n">max_step</span><span class="o">=</span><span class="n">dt_max</span><span class="p">)</span>
    
    <span class="c1"># Push the initial state space to the integrator object</span>
    <span class="n">r</span><span class="o">.</span><span class="n">set_initial_value</span><span class="p">(</span><span class="n">q0</span><span class="p">,</span> <span class="n">t0</span><span class="p">)</span>
    
    <span class="c1"># Make a copy of the q1_local object needed to evaluate the entrainment</span>
    <span class="n">q1_local</span> <span class="o">=</span> <span class="n">deepcopy</span><span class="p">(</span><span class="n">q0_local</span><span class="p">)</span>
    <span class="n">q0_hold</span> <span class="o">=</span> <span class="n">deepcopy</span><span class="p">(</span><span class="n">q1_local</span><span class="p">)</span>
    
    <span class="c1"># Create vectors (using the list data type) to store the solution</span>
    <span class="n">t</span> <span class="o">=</span> <span class="p">[</span><span class="n">t0</span><span class="p">]</span>
    <span class="n">q</span> <span class="o">=</span> <span class="p">[</span><span class="n">q0</span><span class="p">]</span>
    
    <span class="c1"># Integrate a finite number of time steps</span>
    <span class="n">k</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">psteps</span> <span class="o">=</span> <span class="mf">30.</span>
    <span class="n">stop</span> <span class="o">=</span> <span class="kc">False</span>
    <span class="n">neutral_counter</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="n">top_counter</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="k">while</span> <span class="n">r</span><span class="o">.</span><span class="n">successful</span><span class="p">()</span> <span class="ow">and</span> <span class="ow">not</span> <span class="n">stop</span><span class="p">:</span>
        
        <span class="c1"># Print progress to the screen</span>
        <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">remainder</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">float</span><span class="p">(</span><span class="n">k</span><span class="p">),</span> <span class="n">psteps</span><span class="p">)</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
            <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;    Distance:  </span><span class="si">%g</span><span class="s1"> (m), time:  </span><span class="si">%g</span><span class="s1"> (s), k:  </span><span class="si">%d</span><span class="s1">&#39;</span> <span class="o">%</span> \
                <span class="p">(</span><span class="n">q</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="mi">10</span><span class="p">],</span> <span class="n">t</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">],</span> <span class="n">k</span><span class="p">))</span>
        
        <span class="c1"># Perform one step of the integration</span>
        <span class="n">r</span><span class="o">.</span><span class="n">set_f_params</span><span class="p">(</span><span class="n">q0_local</span><span class="p">,</span> <span class="n">q1_local</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">particles</span><span class="p">)</span>
        <span class="n">r</span><span class="o">.</span><span class="n">integrate</span><span class="p">(</span><span class="n">t</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">dt_max</span><span class="p">,</span> <span class="n">step</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
        <span class="n">q1_local</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">r</span><span class="o">.</span><span class="n">t</span><span class="p">,</span> <span class="n">r</span><span class="o">.</span><span class="n">y</span><span class="p">,</span> <span class="n">profile</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">particles</span><span class="p">)</span>
        
        <span class="c1"># Correct the temperature</span>
        <span class="n">r</span> <span class="o">=</span> <span class="n">correct_temperature</span><span class="p">(</span><span class="n">r</span><span class="p">,</span> <span class="n">particles</span><span class="p">)</span>
        
        <span class="c1"># Remove particle solution for particles outside the plume</span>
        <span class="n">r</span> <span class="o">=</span> <span class="n">correct_particle_tracking</span><span class="p">(</span><span class="n">r</span><span class="p">,</span> <span class="n">particles</span><span class="p">)</span>
        
        <span class="c1"># Store the results</span>
        <span class="n">t</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">r</span><span class="o">.</span><span class="n">t</span><span class="p">)</span>
        <span class="n">q</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">r</span><span class="o">.</span><span class="n">y</span><span class="p">)</span>
        
        <span class="c1"># Update the Lagrangian elements for the next time step</span>
        <span class="n">q0_local</span> <span class="o">=</span> <span class="n">q0_hold</span>
        <span class="n">q0_hold</span> <span class="o">=</span> <span class="n">deepcopy</span><span class="p">(</span><span class="n">q1_local</span><span class="p">)</span>
        
        <span class="c1"># Check if the plume has reached a maximum rise height yet</span>
        <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">q0_local</span><span class="o">.</span><span class="n">Jz</span><span class="p">)</span> <span class="o">!=</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">Jz</span><span class="p">):</span>
            <span class="n">top_counter</span> <span class="o">+=</span> <span class="mi">1</span>
        
        <span class="c1"># Check if the plume is at neutral buoyancy in an intrusion layer</span>
        <span class="c1"># (e.g., after the top of the plume)</span>
        <span class="k">if</span> <span class="n">top_counter</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">q0_local</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">rho</span><span class="p">)</span> <span class="o">!=</span> \
                <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">-</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">rho</span><span class="p">):</span>
                <span class="c1"># Update neutral buoyancy level counter</span>
                <span class="n">neutral_counter</span> <span class="o">+=</span> <span class="mi">1</span>
        
        <span class="c1"># Evaluate the stop criteria</span>
        <span class="k">if</span> <span class="n">neutral_counter</span> <span class="o">&gt;=</span> <span class="mi">1</span><span class="p">:</span>
            <span class="c1"># Passed through the second neutral buoyany level</span>
            <span class="n">stop</span> <span class="o">=</span> <span class="kc">True</span>
        <span class="k">if</span> <span class="n">q</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="mi">10</span><span class="p">]</span> <span class="o">/</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">D</span> <span class="o">&gt;</span> <span class="n">sd_max</span><span class="p">:</span>
            <span class="c1"># Progressed desired distance along the plume centerline</span>
            <span class="n">stop</span> <span class="o">=</span> <span class="kc">True</span>
        <span class="k">if</span> <span class="n">k</span> <span class="o">&gt;=</span> <span class="mi">50000</span><span class="p">:</span>
            <span class="c1"># Stop after specified number of iterations; used to protect </span>
            <span class="c1"># against problems with the solution become stuck</span>
            <span class="n">stop</span> <span class="o">=</span> <span class="kc">True</span>
        <span class="k">if</span> <span class="n">q</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="mi">9</span><span class="p">]</span> <span class="o">&lt;=</span> <span class="mf">0.</span><span class="p">:</span>
            <span class="c1"># Reached a location at or above the free surface</span>
            <span class="n">stop</span> <span class="o">=</span> <span class="kc">True</span>
        <span class="k">if</span> <span class="n">q</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="mi">10</span><span class="p">]</span> <span class="o">==</span> <span class="n">q</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">][</span><span class="mi">10</span><span class="p">]:</span>
            <span class="c1"># Progress of motion of the plume has stopped</span>
            <span class="n">stop</span> <span class="o">=</span> <span class="kc">True</span>
        
        <span class="c1"># Update the index counter</span>
        <span class="n">k</span> <span class="o">+=</span> <span class="mi">1</span>
    
    <span class="c1"># Convert solution to numpy arrays</span>
    <span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">t</span><span class="p">)</span>
    <span class="n">q</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">q</span><span class="p">)</span>
    
    <span class="c1"># Show user the final calculated point and return the solution</span>
    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;    Distance:  </span><span class="si">%g</span><span class="s1"> (m), time:  </span><span class="si">%g</span><span class="s1"> (s), k:  </span><span class="si">%d</span><span class="s1">&#39;</span> <span class="o">%</span> \
                <span class="p">(</span><span class="n">q</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span><span class="mi">10</span><span class="p">],</span> <span class="n">t</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">],</span> <span class="n">k</span><span class="p">))</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">q</span><span class="p">)</span></div>


<div class="viewcode-block" id="correct_temperature"><a class="viewcode-back" href="../autodoc/bpm/lmp.correct_temperature.html#lmp.correct_temperature">[docs]</a><span class="k">def</span> <span class="nf">correct_temperature</span><span class="p">(</span><span class="n">r</span><span class="p">,</span> <span class="n">particles</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Make sure the correct temperature is stored in the state space solution</span>
<span class="sd">    </span>
<span class="sd">    When the dispersed phase particles equilibrate to their surrounding</span>
<span class="sd">    temperature, heat transfer is turned off by the methods in </span>
<span class="sd">    `dispersed_phases.Particle`.  This is needed to prevent numerical </span>
<span class="sd">    oscillation as the particle becomes small.  Unfortunately, it is not as</span>
<span class="sd">    easy to make the numerical solution compute the correct result once</span>
<span class="sd">    particle temperature effectively stops being a state space variable since</span>
<span class="sd">    the state space is intrinsic to the ODE solver.  The derivatives function</span>
<span class="sd">    computes the correct heat transfer based on the correct state space, but</span>
<span class="sd">    the state space in the ODE solver remains fixed.</span>
<span class="sd">    </span>
<span class="sd">    Since the solution for heat in the state space of the ODE solver is the</span>
<span class="sd">    wrong value, we have to change the external version of the state space</span>
<span class="sd">    before saving the solution to the current model step.  This follows the</span>
<span class="sd">    same method and reasoning as the similar function in </span>
<span class="sd">    `smp.correct_temperature`.</span>
<span class="sd">    </span>
<span class="sd">    Hence, the purpose of this function is to overwrite the state space </span>
<span class="sd">    solution containing the particle heat that is extrinsic to the ODE solver</span>
<span class="sd">    and which is used to store the state space following each time step.</span>
<span class="sd">    The allows the correct temperature to be stored in the model solution.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    r : `scipy.integrate.ode` object</span>
<span class="sd">        ODE solution containing the current values of the state space in </span>
<span class="sd">        the solver&#39;s extrinsic data.  These values are editable, but an </span>
<span class="sd">        intrinsic version of these data are used when the solver makes </span>
<span class="sd">        calculations; hence, editing this file does not change the state</span>
<span class="sd">        space stored in the actual solver.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `bent_plume_model.Particle` objects containing the dispersed </span>
<span class="sd">        phase local conditions and behavior.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    r : `sciply.integrate.ode` object</span>
<span class="sd">        The updated extrinsic state space with the correct values for heat</span>
<span class="sd">        as were used in the calcualtion.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Find the heat conservation equation in the model state space for the </span>
    <span class="c1"># particles and replace r.y with the correct values.</span>
    <span class="n">idx</span> <span class="o">=</span> <span class="mi">11</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">)):</span>
        <span class="n">idx</span> <span class="o">+=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">nc</span>
        <span class="n">r</span><span class="o">.</span><span class="n">y</span><span class="p">[</span><span class="n">idx</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">m</span><span class="p">)</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nbe</span> <span class="o">*</span> \
                       <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">cp</span> <span class="o">*</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">T</span>
        <span class="c1"># Advance for heat, time, and position</span>
        <span class="n">idx</span> <span class="o">+=</span> <span class="mi">1</span> <span class="o">+</span> <span class="mi">1</span> <span class="o">+</span> <span class="mi">3</span>
    
    <span class="c1"># Return the corrected solution</span>
    <span class="k">return</span> <span class="n">r</span></div>


<div class="viewcode-block" id="correct_particle_tracking"><a class="viewcode-back" href="../autodoc/bpm/lmp.correct_particle_tracking.html#lmp.correct_particle_tracking">[docs]</a><span class="k">def</span> <span class="nf">correct_particle_tracking</span><span class="p">(</span><span class="n">r</span><span class="p">,</span> <span class="n">particles</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Remove the particle tracking solution after particles exit plume</span>
<span class="sd">    </span>
<span class="sd">    Even though the particle tracking stops as needed once the particles</span>
<span class="sd">    leave the plume, the post processing algorithm has now way to know if a</span>
<span class="sd">    given state space solution is before or after particle tracking has</span>
<span class="sd">    stopped. This function simply replaces the particle position after</span>
<span class="sd">    integration has stopped (e.g., after the particles leave the plume) with</span>
<span class="sd">    NaN so that the post-processor always knows whether the solution in the</span>
<span class="sd">    state space is valid or not.  This is necessary since the solution for </span>
<span class="sd">    particle position is in local plume coordinates (l,n,m); hence, it is not</span>
<span class="sd">    possible to know the (x,y,z) position unless the correct local plume </span>
<span class="sd">    element is known.  This function makes sure that every valid (l,n,m) is </span>
<span class="sd">    stored with the corresponding element.  </span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    r : `scipy.integrate.ode` object</span>
<span class="sd">        ODE solution containing the current values of the state space in </span>
<span class="sd">        the solver&#39;s extrinsic data.  These values are editable, but an </span>
<span class="sd">        intrinsic version of these data are used when the solver makes </span>
<span class="sd">        calculations; hence, editing this file does not change the state</span>
<span class="sd">        space stored in the actual solver.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `bent_plume_model.Particle` objects containing the dispersed </span>
<span class="sd">        phase local conditions and behavior.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    r : `sciply.integrate.ode` object</span>
<span class="sd">        The updated extrinsic state space with the correct values for heat</span>
<span class="sd">        as were used in the calcualtion.</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Skip through the single-phase state space</span>
    <span class="n">idx</span> <span class="o">=</span> <span class="mi">11</span>

    <span class="c1"># Check each particle to determine whether they are inside or outside</span>
    <span class="c1"># the plume</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">)):</span>
        <span class="k">if</span> <span class="ow">not</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">integrate</span><span class="p">:</span>
            <span class="c1"># Skip the masses, temperature, and time</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">nc</span> <span class="o">+</span> <span class="mi">2</span>
            
            <span class="c1"># Particle is outside the plume; replace the coordinates with </span>
            <span class="c1"># np.nan</span>
            <span class="n">r</span><span class="o">.</span><span class="n">y</span><span class="p">[</span><span class="n">idx</span><span class="p">:</span><span class="n">idx</span><span class="o">+</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">nan</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="mi">3</span>
            
        <span class="k">else</span><span class="p">:</span>
            <span class="c1"># Skip the masses, temperature, time, and coordinates</span>
            <span class="n">idx</span> <span class="o">+=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">particle</span><span class="o">.</span><span class="n">nc</span> <span class="o">+</span> <span class="mi">5</span>
        
        <span class="c1"># Check if the integration should stop</span>
        <span class="k">if</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">p_fac</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
            <span class="c1"># Stop tracking the particle inside the plume</span>
            <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">integrate</span> <span class="o">=</span> <span class="kc">False</span>
            
            <span class="c1"># Store the properties at the exit point</span>
            <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">te</span> <span class="o">=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">t</span>
            <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">xe</span> <span class="o">=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">x</span>
            <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">ye</span> <span class="o">=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">y</span>
            <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">ze</span> <span class="o">=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">z</span>
            <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">me</span> <span class="o">=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">m</span>
            <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">Te</span> <span class="o">=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">T</span>
    
    <span class="c1"># Return the corrected solution</span>
    <span class="k">return</span> <span class="n">r</span></div>


<div class="viewcode-block" id="entrainment"><a class="viewcode-back" href="../autodoc/bpm/lmp.entrainment.html#lmp.entrainment">[docs]</a><span class="k">def</span> <span class="nf">entrainment</span><span class="p">(</span><span class="n">q0_local</span><span class="p">,</span> <span class="n">q1_local</span><span class="p">,</span> <span class="n">p</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the total shear and forced entrainment at one time step</span>
<span class="sd">    </span>
<span class="sd">    Computes the local entrainment (kg/s) as a combination of shear</span>
<span class="sd">    entrainment and forced entrainment for a local Lagrangian element. This</span>
<span class="sd">    function follows the approach in Lee and Cheung (1990) to compute both</span>
<span class="sd">    types of entrainment, but uses the formulation in Jirka (2004) for the</span>
<span class="sd">    shear entrainment term. Like Lee and Cheung (1990), it uses the maximum</span>
<span class="sd">    entrainment hypothesis: entrainment = max (shear, forced), with the</span>
<span class="sd">    exception that a pure coflowing momentum jet has entrainment = shear +</span>
<span class="sd">    forced. This function also makes one correction that in pure coflow</span>
<span class="sd">    the forced entrainment should be computed by integrating around the entire</span>
<span class="sd">    jet, and not just the half of the jet exposed to the current.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    q0_local : `bent_plume_model.LagElement`</span>
<span class="sd">        Object containing the numerical solution at the previous time step</span>
<span class="sd">    q1_local : `bent_plume_model.LagElement`</span>
<span class="sd">        Object containing the numerical solution at the current time step</span>
<span class="sd">    p : `ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the bent</span>
<span class="sd">        plume model.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    md : float</span>
<span class="sd">        Total entrainment (kg/s)</span>
<span class="sd">    </span>
<span class="sd">    Notes</span>
<span class="sd">    -----</span>
<span class="sd">    The entrainment computed here is already integrated over the current </span>
<span class="sd">    Lagrangian element surface area.  Hence, the result is (kg/s) into the </span>
<span class="sd">    element.  </span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Find the magnitude and direction of the velocity vector in q1_local</span>
    <span class="n">Ua</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">ua</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">va</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">wa</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
    <span class="n">Phi_a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arctan2</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">wa</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">ua</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">va</span><span class="o">**</span><span class="mi">2</span><span class="p">))</span>
    <span class="n">Theta_a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arctan2</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">va</span><span class="p">,</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">ua</span><span class="p">)</span>
    
    <span class="c1"># Get the component of the ambient current along the plume centerline</span>
    <span class="n">Us</span> <span class="o">=</span> <span class="n">Ua</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">phi</span> <span class="o">-</span> <span class="n">Phi_a</span><span class="p">)</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">-</span> <span class="n">Theta_a</span><span class="p">)</span>
    
    <span class="c1"># Get the sines and cosines of the new angles</span>
    <span class="n">sin_t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">-</span> <span class="n">Theta_a</span><span class="p">)</span>
    <span class="n">sin_p</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">phi</span> <span class="o">-</span> <span class="n">Phi_a</span><span class="p">)</span>
    <span class="n">cos_t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">-</span> <span class="n">Theta_a</span><span class="p">)</span>
    <span class="n">cos_p</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">phi</span> <span class="o">-</span> <span class="n">Phi_a</span><span class="p">)</span>
    <span class="n">cos_t0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">q0_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">-</span> <span class="n">Theta_a</span><span class="p">)</span>
    <span class="n">cos_p0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">q0_local</span><span class="o">.</span><span class="n">phi</span> <span class="o">-</span> <span class="n">Phi_a</span><span class="p">)</span>
    
    <span class="c1"># Get the shear entrainment coefficient for the top-hat model.  In this</span>
    <span class="c1"># equation, phi has to be with reference to the gravity vector; hence, </span>
    <span class="c1"># we pass phi for the fixed coordinate system, but theta has to be the</span>
    <span class="c1"># angle from the crossflow direction, so be pass theta - theta_a.</span>
    <span class="n">alpha_s</span> <span class="o">=</span> <span class="n">dispersed_phases</span><span class="o">.</span><span class="n">shear_entrainment</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">V</span><span class="p">,</span> <span class="n">Us</span><span class="p">,</span> 
        <span class="n">q1_local</span><span class="o">.</span><span class="n">rho</span><span class="p">,</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">rho_a</span><span class="p">,</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">b</span><span class="p">,</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_p</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
    
    <span class="c1"># Total shear entrainment (kg/s)</span>
    <span class="n">md_s</span> <span class="o">=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">V</span> <span class="o">-</span> <span class="n">Us</span><span class="p">)</span> <span class="o">*</span> <span class="n">alpha_s</span> <span class="o">*</span> <span class="p">(</span> <span class="mf">2.</span> 
           <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">h</span><span class="p">)</span>
         
    <span class="c1"># Compute the projected area entrainment terms...first, the crossflow </span>
    <span class="c1"># projected onto the plume centerline</span>
    <span class="n">a1</span> <span class="o">=</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">sin_p</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">sin_t</span><span class="o">**</span><span class="mi">2</span> <span class="o">-</span> <span class="n">sin_p</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> 
         <span class="n">sin_t</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">h</span>
    <span class="k">if</span> <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">s</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">s</span><span class="p">)</span> <span class="o">/</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">b</span> <span class="o">&lt;=</span> <span class="mf">1.e-3</span><span class="p">:</span>
        <span class="c1"># The plume is not progressing along the centerline; assume the</span>
        <span class="c1"># expansion and curvature corrections are small since delta s / b is</span>
        <span class="c1"># very small.</span>
        <span class="n">a2</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="n">a3</span> <span class="o">=</span> <span class="mf">0.</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># Second, correction for plume expansion</span>
        <span class="n">a2</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">b</span> <span class="o">*</span> <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">b</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">b</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">s</span> <span class="o">-</span> 
             <span class="n">q0_local</span><span class="o">.</span><span class="n">s</span><span class="p">)</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">h</span> <span class="o">*</span> <span class="n">cos_p</span> <span class="o">*</span> <span class="n">cos_t</span>
        <span class="c1"># Third, correction for plume curvature</span>
        <span class="n">a3</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">b</span><span class="o">**</span><span class="mi">2</span> <span class="o">/</span> <span class="mf">2.</span> <span class="o">*</span> <span class="p">(</span><span class="n">cos_p</span> <span class="o">*</span> <span class="n">cos_t</span> <span class="o">-</span> <span class="n">cos_p0</span> <span class="o">*</span> 
             <span class="n">cos_t0</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">s</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">s</span><span class="p">)</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">h</span>
    
    <span class="c1"># Get the total projected area for the forced entraiment</span>
    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">sin_t</span><span class="p">)</span> <span class="o">&lt;=</span> <span class="mf">1.e-9</span> <span class="ow">and</span> <span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">sin_p</span><span class="p">)</span> <span class="o">&lt;=</span> <span class="mf">1.e-9</span><span class="p">:</span>
        <span class="c1"># Jet is in co-flow, shear entrainment model takes care of this case</span>
        <span class="c1"># by itself</span>
        <span class="n">A</span> <span class="o">=</span> <span class="mf">0.</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">A</span> <span class="o">=</span> <span class="n">a1</span> <span class="o">+</span> <span class="n">a2</span> <span class="o">+</span> <span class="n">a3</span>
    
    <span class="c1"># Total forced entrainment (kg/s)</span>
    <span class="n">md_f</span> <span class="o">=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">*</span> <span class="n">Ua</span> <span class="o">*</span> <span class="n">A</span>
    
    <span class="c1"># Obtain the total entrainment using the maximum hypothesis from Lee and </span>
    <span class="c1"># Cheung (1990)</span>
    <span class="k">if</span> <span class="n">md_s</span> <span class="o">&gt;</span> <span class="n">md_f</span><span class="p">:</span>
        <span class="n">md</span> <span class="o">=</span> <span class="n">md_s</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">md</span> <span class="o">=</span> <span class="n">md_f</span>
    
    <span class="c1"># Return the entrainment rate</span>
    <span class="k">return</span> <span class="n">md</span></div>


<div class="viewcode-block" id="track_particles"><a class="viewcode-back" href="../autodoc/bpm/lmp.track_particles.html#lmp.track_particles">[docs]</a><span class="k">def</span> <span class="nf">track_particles</span><span class="p">(</span><span class="n">q0_local</span><span class="p">,</span> <span class="n">q1_local</span><span class="p">,</span> <span class="n">md</span><span class="p">,</span> <span class="n">particles</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the forcing variables needed to track particles</span>
<span class="sd">    </span>
<span class="sd">    The independent variable in the Lagrangian plume model is t for advection</span>
<span class="sd">    of the continuous phase.  Because the particles slip through the fluid, </span>
<span class="sd">    their advection speed is different; hence, all particle equations have </span>
<span class="sd">    a different independent variable for their time, tp.  Also, particle </span>
<span class="sd">    motion is computed in the local plume coordinates space (l,n,m); thus,</span>
<span class="sd">    the vertical slip velocity needs to be transformed to the local plume</span>
<span class="sd">    coordinate system.  Finally, the entrainment velocity pointing toward</span>
<span class="sd">    the plume centerline needs to be evaluated, based on the distance the </span>
<span class="sd">    particle is from the plume centerline.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    q0_local : `bent_plume_model.LagElement`</span>
<span class="sd">        Object containing the numerical solution at the previous time step</span>
<span class="sd">    q1_local : `bent_plume_model.LagElement`</span>
<span class="sd">        Object containing the numerical solution at the current time step</span>
<span class="sd">    md : float</span>
<span class="sd">        Total entrainment into the Lagrangian element (kg)</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `bent_plume_model.Particle` objects containing the dispersed </span>
<span class="sd">        phase local conditions and behavior.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    fe : float</span>
<span class="sd">        Entrainment frequency (1/s)</span>
<span class="sd">    up : ndarray</span>
<span class="sd">        Slip velocity for each particle projected on the local plume </span>
<span class="sd">        coordinate system (l,n,m) (m/s).  Each row is for a different </span>
<span class="sd">        particle and the columns are for the velocity in the (l,n,m)</span>
<span class="sd">        direction.</span>
<span class="sd">    dtp_dt : ndarray</span>
<span class="sd">        Differential particle transport time to continuous phase transport</span>
<span class="sd">        time evaluated from the previous time step to the current time </span>
<span class="sd">        step (--).</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Compute the entrainment frequency</span>
    <span class="n">fe</span> <span class="o">=</span> <span class="n">md</span> <span class="o">/</span> <span class="p">(</span><span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">rho_a</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">b</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> 
         <span class="n">q1_local</span><span class="o">.</span><span class="n">h</span><span class="p">)</span>
    
    <span class="c1"># Get the rotation matrix to the local coordinate system (l,n,m)</span>
    <span class="n">ds</span> <span class="o">=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">s</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">s</span>
    <span class="n">A</span> <span class="o">=</span> <span class="n">local_coords</span><span class="p">(</span><span class="n">q0_local</span><span class="p">,</span> <span class="n">q1_local</span><span class="p">,</span> <span class="n">ds</span><span class="p">)</span>
    
    <span class="c1"># Get the velocity of the current plume slice</span>
    <span class="n">V</span> <span class="o">=</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">V</span>
    
    <span class="c1"># Compute particle properties</span>
    <span class="n">up</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">),</span><span class="mi">3</span><span class="p">))</span>
    <span class="n">dtp_dt</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">))</span>
    
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">)):</span>
        <span class="c1"># Transform the slip velocity from Cartesian coordinate to the </span>
        <span class="c1"># local plume coordinate system (l,n,m)</span>
        <span class="n">up</span><span class="p">[</span><span class="n">i</span><span class="p">,:]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mf">0.</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="o">-</span><span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">us</span><span class="p">]))</span>
        
        <span class="c1"># Get the distance along the particle path</span>
        <span class="n">dsp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">q1_local</span><span class="o">.</span><span class="n">x_p</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">x_p</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">0</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> 
                      <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">x_p</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">x_p</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">1</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> 
                      <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">x_p</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">x_p</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">2</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
        
        <span class="c1"># Get the total velocity of the particle</span>
        <span class="n">Vp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">up</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> 
                     <span class="p">(</span><span class="n">up</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">fe</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">X_p</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">1</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> 
                     <span class="p">(</span><span class="n">up</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">-</span> <span class="n">fe</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">X_p</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="mi">2</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
        
        <span class="c1"># Compute the particle time correction dtp/dt</span>
        <span class="k">if</span> <span class="n">Vp</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="k">elif</span> <span class="n">ds</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">dtp_dt</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">V</span> <span class="o">/</span> <span class="n">Vp</span> <span class="o">*</span> <span class="n">dsp</span> <span class="o">/</span> <span class="n">ds</span>
    
    <span class="c1"># Return the particle tracking variables</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">fe</span><span class="p">,</span> <span class="n">up</span><span class="p">,</span> <span class="n">dtp_dt</span><span class="p">)</span></div>


<div class="viewcode-block" id="local_coords"><a class="viewcode-back" href="../autodoc/bpm/lmp.local_coords.html#lmp.local_coords">[docs]</a><span class="k">def</span> <span class="nf">local_coords</span><span class="p">(</span><span class="n">q0_local</span><span class="p">,</span> <span class="n">q1_local</span><span class="p">,</span> <span class="n">ds</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the rotation matrix from (x, y, z) to (l, n, m)</span>
<span class="sd">    </span>
<span class="sd">    Computes the rotation matrix from the local Cartesian coordinate system</span>
<span class="sd">    (x - xi, y - yi, z - zi), where (xi, yi, zi) is the current location of</span>
<span class="sd">    the Lagrangian plume element, to the system tangent to the current plume</span>
<span class="sd">    trajectory (l, n, m); l is oriented tangent to the plume centerline, </span>
<span class="sd">    n is orthogonal to l and along the line from the local radius of </span>
<span class="sd">    curvature, and m is orthogonal to n and l.  The transformation is </span>
<span class="sd">    provided in Lee and Chueng (1990).  This method makes a small adaptation</span>
<span class="sd">    allowing for infinite radius of curvature (plume propagating along a </span>
<span class="sd">    straight line).  </span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    q0_local : `bent_plume_model.LagElement`</span>
<span class="sd">        Object containing the numerical solution at the previous time step</span>
<span class="sd">    q1_local : `bent_plume_model.LagElement`</span>
<span class="sd">        Object containing the numerical solution at the current time step</span>
<span class="sd">    ds : float</span>
<span class="sd">        Segment length along the centerline between the solutions contained </span>
<span class="sd">        in `q0_local` and `q1_local`.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    A : ndarray</span>
<span class="sd">        Rotation matrix from (x, y, z)&#39; to (l, n, m).  The inverse of this</span>
<span class="sd">        matrix will convert back from (l, n, m) to (x, y, z)&#39;.</span>
<span class="sd">    </span>
<span class="sd">    See Also</span>
<span class="sd">    --------</span>
<span class="sd">    bent_plume_model.Particle.track</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the rate of angular rotation for the centerline</span>
    <span class="k">if</span> <span class="n">ds</span> <span class="o">&lt;</span> <span class="mf">1.e-12</span><span class="p">:</span>
        <span class="n">phi_d</span> <span class="o">=</span> <span class="mf">0.</span>
        <span class="n">theta_d</span> <span class="o">=</span> <span class="mf">0.</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">phi_d</span> <span class="o">=</span> <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">phi</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">phi</span><span class="p">)</span> <span class="o">/</span> <span class="n">ds</span>
        <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">theta</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">:</span>
            <span class="c1"># Angles are close to 0</span>
            <span class="k">if</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">&gt;</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">:</span>
                <span class="n">theta_d</span> <span class="o">=</span> <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">-</span> <span class="p">(</span><span class="n">q0_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">+</span> <span class="mf">2.</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span><span class="p">))</span> <span class="o">/</span> <span class="n">ds</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">theta_d</span> <span class="o">=</span> <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">+</span> <span class="mf">2.</span><span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">theta</span><span class="p">)</span> <span class="o">/</span> <span class="n">ds</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">theta_d</span> <span class="o">=</span> <span class="p">(</span><span class="n">q1_local</span><span class="o">.</span><span class="n">theta</span> <span class="o">-</span> <span class="n">q0_local</span><span class="o">.</span><span class="n">theta</span><span class="p">)</span> <span class="o">/</span> <span class="n">ds</span>
    
    <span class="c1"># Get the value of 1 / R = r</span>
    <span class="n">r</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">phi_d</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">theta_d</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
    
    <span class="c1"># Compute the rotation matrix between (i,j,k) and (l,n,m)</span>
    <span class="k">if</span> <span class="n">r</span> <span class="o">&lt;</span> <span class="mf">1.e-8</span><span class="p">:</span>
        <span class="c1"># Trajectory is straight and R is infinite</span>
        <span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span>
            <span class="p">[</span><span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_t</span><span class="p">,</span> 
             <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_t</span><span class="p">,</span> 
             <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_p</span><span class="p">],</span> 
            <span class="p">[</span><span class="n">q1_local</span><span class="o">.</span><span class="n">cos_t</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_p</span><span class="p">,</span>
             <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_t</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_p</span><span class="p">,</span>
             <span class="o">-</span><span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span><span class="p">],</span> 
            <span class="p">[</span><span class="n">q1_local</span><span class="o">.</span><span class="n">sin_t</span><span class="p">,</span>
             <span class="o">-</span><span class="n">q1_local</span><span class="o">.</span><span class="n">cos_t</span><span class="p">,</span>
             <span class="mf">0.</span><span class="p">]])</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># Trajectory is curving, and R is finite</span>
        <span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span>
            <span class="p">[</span><span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_t</span><span class="p">,</span> 
             <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_t</span><span class="p">,</span> 
             <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_p</span><span class="p">],</span> 
            <span class="p">[</span><span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_t</span> <span class="o">*</span> <span class="n">theta_d</span> <span class="o">+</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_t</span> <span class="o">*</span> 
                <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_p</span> <span class="o">*</span> <span class="n">phi_d</span><span class="p">,</span> 
             <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_t</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_p</span> <span class="o">*</span> <span class="n">phi_d</span> <span class="o">-</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span> <span class="o">*</span> 
                <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_t</span> <span class="o">*</span> <span class="n">theta_d</span><span class="p">,</span> 
             <span class="o">-</span><span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span> <span class="o">*</span> <span class="n">phi_d</span><span class="p">]</span> <span class="o">/</span> <span class="n">r</span><span class="p">,</span> 
            <span class="p">[</span><span class="n">q1_local</span><span class="o">.</span><span class="n">sin_t</span> <span class="o">*</span> <span class="n">phi_d</span> <span class="o">-</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_p</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span> <span class="o">*</span> 
                <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_t</span> <span class="o">*</span> <span class="n">theta_d</span><span class="p">,</span>
             <span class="o">-</span><span class="n">q1_local</span><span class="o">.</span><span class="n">cos_t</span> <span class="o">*</span> <span class="n">phi_d</span> <span class="o">-</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span> <span class="o">*</span> <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_t</span> <span class="o">*</span> 
             <span class="n">q1_local</span><span class="o">.</span><span class="n">sin_p</span> <span class="o">*</span> <span class="n">theta_d</span><span class="p">,</span> 
             <span class="n">q1_local</span><span class="o">.</span><span class="n">cos_p</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">theta_d</span><span class="p">]</span> <span class="o">/</span> <span class="n">r</span><span class="p">])</span>
    
    <span class="c1"># Return the rotation matrix</span>
    <span class="k">return</span> <span class="n">A</span></div>


<span class="c1"># ----------------------------------------------------------------------------</span>
<span class="c1"># Functions to compute the initial conditions for the first model element</span>
<span class="c1"># ----------------------------------------------------------------------------</span>

<div class="viewcode-block" id="main_ic"><a class="viewcode-back" href="../autodoc/bpm/lmp.main_ic.html#lmp.main_ic">[docs]</a><span class="k">def</span> <span class="nf">main_ic</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">X</span><span class="p">,</span> <span class="n">D</span><span class="p">,</span> <span class="n">Vj</span><span class="p">,</span> <span class="n">phi_0</span><span class="p">,</span> <span class="n">theta_0</span><span class="p">,</span> <span class="n">Sj</span><span class="p">,</span> <span class="n">Tj</span><span class="p">,</span> <span class="n">cj</span><span class="p">,</span> 
            <span class="n">tracers</span><span class="p">,</span> <span class="n">p</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the initial conditions for the Lagrangian plume state space</span>
<span class="sd">    </span>
<span class="sd">    Compute the initial conditions at the release location for a Lagrangian</span>
<span class="sd">    plume element.  This can either be a pure single-phase plume, a pure</span>
<span class="sd">    multiphase plume, or a mixed release of multiphase and continuous phase</span>
<span class="sd">    fluid.</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the single bubble model simulation.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `bent_plume_model.Particle` objects containing the dispersed </span>
<span class="sd">        phase local conditions and behavior.</span>
<span class="sd">    X : ndarray</span>
<span class="sd">        Release location (x, y, z) in (m)</span>
<span class="sd">    D : float</span>
<span class="sd">        Diameter for the equivalent circular cross-section of the release </span>
<span class="sd">        (m)</span>
<span class="sd">    Vj : float</span>
<span class="sd">        Scalar value of the magnitude of the discharge velocity for </span>
<span class="sd">        continuous phase fluid in the discharge.  This variable should be </span>
<span class="sd">        0 or None for a pure multiphase discharge.</span>
<span class="sd">    phi_0 : float</span>
<span class="sd">        Vertical angle from the horizontal for the discharge orientation </span>
<span class="sd">        (rad in range +/- pi/2)</span>
<span class="sd">    theta_0 : float</span>
<span class="sd">        Horizontal angle from the x-axis for the discharge orientation.  </span>
<span class="sd">        The x-axis is taken in the direction of the ambient current.  </span>
<span class="sd">        (rad in range 0 to 2 pi)</span>
<span class="sd">    Sj : float</span>
<span class="sd">        Salinity of the continuous phase fluid in the discharge (psu)</span>
<span class="sd">    Tj : float</span>
<span class="sd">        Temperature of the continuous phase fluid in the discharge (T)</span>
<span class="sd">    cj : ndarray</span>
<span class="sd">        Concentration of passive tracers in the discharge (user-defined)</span>
<span class="sd">    tracers : string list</span>
<span class="sd">        List of passive tracers in the discharge.  These can be chemicals </span>
<span class="sd">        present in the ambient `profile` data, and if so, entrainment of </span>
<span class="sd">        these chemicals will change the concentrations computed for these </span>
<span class="sd">        tracers.  However, none of these concentrations are used in the </span>
<span class="sd">        dissolution of the dispersed phase.  Hence, `tracers` should not </span>
<span class="sd">        contain any chemicals present in the dispersed phase particles.</span>
<span class="sd">    p : `stratified_plume_model.ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    t : float</span>
<span class="sd">        Initial time for the simulation (s)</span>
<span class="sd">    q : ndarray</span>
<span class="sd">        Initial value of the plume state space</span>
<span class="sd">    chem_names : str list</span>
<span class="sd">        List of the chemicals in the dispersed phase composition that are </span>
<span class="sd">        undergoing dissolution</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the initial volume flux</span>
    <span class="k">if</span> <span class="n">Vj</span> <span class="ow">is</span> <span class="kc">None</span> <span class="ow">or</span> <span class="n">Vj</span> <span class="o">==</span> <span class="mf">0.</span><span class="p">:</span>
        <span class="c1"># This is a pure multiphase plume.  Estimate the initial conditions </span>
        <span class="c1"># using Wuest et al. 1992.</span>
        <span class="n">Q</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">X</span><span class="p">,</span> <span class="n">Tj</span><span class="p">,</span> <span class="n">Sj</span><span class="p">,</span> <span class="n">Pj</span><span class="p">,</span> <span class="n">rho_j</span> <span class="o">=</span> <span class="n">dispersed_phases</span><span class="o">.</span><span class="n">zfe_volume_flux</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> 
            <span class="n">particles</span><span class="p">,</span> <span class="n">p</span><span class="p">,</span> <span class="n">X</span><span class="p">,</span> <span class="n">D</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span>
    
    <span class="k">else</span><span class="p">:</span>
        <span class="c1"># The discharge contains continuous phase fluid.  Get the flow rate</span>
        <span class="c1"># from the discharge conditions.</span>
        <span class="n">Q</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">X</span><span class="p">,</span> <span class="n">Tj</span><span class="p">,</span> <span class="n">Sj</span><span class="p">,</span> <span class="n">Pj</span><span class="p">,</span> <span class="n">rho_j</span> <span class="o">=</span> <span class="n">zfe_volume_flux</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> <span class="n">X</span><span class="p">,</span> <span class="n">D</span><span class="o">/</span><span class="mf">2.</span><span class="p">,</span> <span class="n">Vj</span><span class="p">,</span> 
                                     <span class="n">Sj</span><span class="p">,</span> <span class="n">Tj</span><span class="p">)</span>
    
    <span class="c1"># Get the names of the chemicals to track</span>
    <span class="n">chem_names</span> <span class="o">=</span> <span class="n">dispersed_phases</span><span class="o">.</span><span class="n">get_chem_names</span><span class="p">(</span><span class="n">particles</span><span class="p">)</span>
    
    <span class="c1"># Build the initial state space with these initial values</span>
    <span class="n">t</span><span class="p">,</span> <span class="n">q</span> <span class="o">=</span> <span class="n">bent_plume_ic</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">Q</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">D</span><span class="p">,</span> <span class="n">X</span><span class="p">,</span> <span class="n">phi_0</span><span class="p">,</span> <span class="n">theta_0</span><span class="p">,</span> 
             <span class="n">Tj</span><span class="p">,</span> <span class="n">Sj</span><span class="p">,</span> <span class="n">Pj</span><span class="p">,</span> <span class="n">rho_j</span><span class="p">,</span> <span class="n">cj</span><span class="p">,</span> <span class="n">chem_names</span><span class="p">,</span> <span class="n">tracers</span><span class="p">,</span> <span class="n">p</span><span class="p">)</span>
    
    <span class="c1"># Return the initial state space</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">q</span><span class="p">,</span> <span class="n">chem_names</span><span class="p">)</span></div>


<div class="viewcode-block" id="bent_plume_ic"><a class="viewcode-back" href="../autodoc/bpm/lmp.bent_plume_ic.html#lmp.bent_plume_ic">[docs]</a><span class="k">def</span> <span class="nf">bent_plume_ic</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> <span class="n">particles</span><span class="p">,</span> <span class="n">Qj</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">D</span><span class="p">,</span> <span class="n">X</span><span class="p">,</span> <span class="n">phi_0</span><span class="p">,</span> <span class="n">theta_0</span><span class="p">,</span> <span class="n">Tj</span><span class="p">,</span> <span class="n">Sj</span><span class="p">,</span> 
                  <span class="n">Pj</span><span class="p">,</span> <span class="n">rho_j</span><span class="p">,</span> <span class="n">cj</span><span class="p">,</span> <span class="n">chem_names</span><span class="p">,</span> <span class="n">tracers</span><span class="p">,</span> <span class="n">p</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Build the Lagragian plume state space given the initial conditions</span>
<span class="sd">    </span>
<span class="sd">    Constructs the initial state space for a Lagrangian plume element from </span>
<span class="sd">    the initial values for the base plume variables (e.g., Q, J, u, S, T, </span>
<span class="sd">    etc.).</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the single bubble model simulation.</span>
<span class="sd">    particles : list of `Particle` objects</span>
<span class="sd">        List of `bent_plume_model.Particle` objects containing the dispersed </span>
<span class="sd">        phase local conditions and behavior.</span>
<span class="sd">    Qj : Volume flux of continuous phase fluid at the discharge (m^3/s)</span>
<span class="sd">    A : Cross-sectional area of the discharge (M^2)</span>
<span class="sd">    D : float</span>
<span class="sd">        Diameter for the equivalent circular cross-section of the release </span>
<span class="sd">        (m)</span>
<span class="sd">    X : ndarray</span>
<span class="sd">        Release location (x, y, z) in (m)</span>
<span class="sd">    phi_0 : float</span>
<span class="sd">        Vertical angle from the horizontal for the discharge orientation </span>
<span class="sd">        (rad in range +/- pi/2)</span>
<span class="sd">    theta_0 : float</span>
<span class="sd">        Horizontal angle from the x-axis for the discharge orientation.  </span>
<span class="sd">        The x-axis is taken in the direction of the ambient current.  </span>
<span class="sd">        (rad in range 0 to 2 pi)</span>
<span class="sd">    Tj : float</span>
<span class="sd">        Temperature of the continuous phase fluid in the discharge (T)</span>
<span class="sd">    Sj : float</span>
<span class="sd">        Salinity of the continuous phase fluid in the discharge (psu)</span>
<span class="sd">    Pj : float</span>
<span class="sd">        Pressure at the discharge (Pa)</span>
<span class="sd">    rho_j : float</span>
<span class="sd">        Density of the continous phase fluid in the discharge (kg/m^3)</span>
<span class="sd">    cj : ndarray</span>
<span class="sd">        Concentration of passive tracers in the discharge (user-defined)</span>
<span class="sd">    chem_names : string list</span>
<span class="sd">        List of chemical parameters to track for the dissolution.  Only the </span>
<span class="sd">        parameters in this list will be used to set background concentration</span>
<span class="sd">        for the dissolution, and the concentrations of these parameters are </span>
<span class="sd">        computed separately from those listed in `tracers` or inputed from</span>
<span class="sd">        the discharge through `cj`.</span>
<span class="sd">    tracers : string list</span>
<span class="sd">        List of passive tracers in the discharge.  These can be chemicals </span>
<span class="sd">        present in the ambient `profile` data, and if so, entrainment of </span>
<span class="sd">        these chemicals will change the concentrations computed for these </span>
<span class="sd">        tracers.  However, none of these concentrations are used in the </span>
<span class="sd">        dissolution of the dispersed phase.  Hence, `tracers` should not </span>
<span class="sd">        contain any chemicals present in the dispersed phase particles.</span>
<span class="sd">    p : `stratified_plume_model.ModelParams` object</span>
<span class="sd">        Object containing the fixed model parameters for the stratified </span>
<span class="sd">        plume model.</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    t : float</span>
<span class="sd">        Initial time for the simulation (s)</span>
<span class="sd">    q : ndarray</span>
<span class="sd">        Initial value of the plume state space</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    
    <span class="c1"># Set the dimensions of the initial Lagrangian plume element.</span>
    <span class="n">b</span> <span class="o">=</span> <span class="n">D</span> <span class="o">/</span> <span class="mf">2.</span>
    <span class="n">h</span> <span class="o">=</span> <span class="n">D</span> <span class="o">/</span> <span class="mf">5.</span>
    
    <span class="c1"># Measure the arc length along the plume</span>
    <span class="n">s0</span> <span class="o">=</span> <span class="mf">0.</span>
    
    <span class="c1"># The total discharge volume flux is the jet discharge since we assume</span>
    <span class="c1"># the void fraction of gas is negligible</span>
    <span class="n">Q</span> <span class="o">=</span> <span class="n">Qj</span>
    
    <span class="c1"># Determine the time to fill the initial Lagrangian element</span>
    <span class="n">dt</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">b</span><span class="o">**</span><span class="mi">2</span> <span class="o">*</span> <span class="n">h</span> <span class="o">/</span> <span class="n">Q</span>
    
    <span class="c1"># Compute the mass of jet discharge in the initial Lagrangian element</span>
    <span class="n">Mj</span> <span class="o">=</span> <span class="n">Qj</span> <span class="o">*</span> <span class="n">dt</span> <span class="o">*</span> <span class="n">rho_j</span>
    
    <span class="c1"># Evaluate the mass of particles in the intial Lagrangian element.  Since</span>
    <span class="c1"># particles are tracked by number and mass per particle, we need to know</span>
    <span class="c1"># how many particles enter the Lagrangian element.  This should be the </span>
    <span class="c1"># number flux in #/s time the fill time for the Lagrangian element, dt.</span>
    <span class="c1"># Store this value in the `Particle` objects for use throughout the model.</span>
    <span class="n">nbe</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">))</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">particles</span><span class="p">)):</span>
        <span class="n">nbe</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nb0</span> <span class="o">*</span> <span class="n">dt</span>
        <span class="n">particles</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">.</span><span class="n">nbe</span> <span class="o">=</span> <span class="n">nbe</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
    
    <span class="c1"># Get the velocity in the component directions</span>
    <span class="n">Uj</span> <span class="o">=</span> <span class="n">flux_to_velocity</span><span class="p">(</span><span class="n">Qj</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">phi_0</span><span class="p">,</span> <span class="n">theta_0</span><span class="p">)</span>
    
    <span class="c1"># Compute the magnitude of the exit velocity </span>
    <span class="n">V</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">Uj</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">Uj</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="n">Uj</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
    
    <span class="c1"># Build the continuous-phase portion of the model state space vector</span>
    <span class="n">t</span> <span class="o">=</span> <span class="mf">0.</span>
    <span class="n">q</span> <span class="o">=</span> <span class="p">[</span><span class="n">Mj</span><span class="p">,</span> <span class="n">Mj</span> <span class="o">*</span> <span class="n">Sj</span><span class="p">,</span> <span class="n">Mj</span> <span class="o">*</span> <span class="n">seawater</span><span class="o">.</span><span class="n">cp</span><span class="p">()</span> <span class="o">*</span> <span class="n">Tj</span><span class="p">,</span> <span class="n">Mj</span> <span class="o">*</span> <span class="n">Uj</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">Mj</span> <span class="o">*</span> <span class="n">Uj</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> 
          <span class="n">Mj</span> <span class="o">*</span> <span class="n">Uj</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="n">h</span> <span class="o">/</span> <span class="n">V</span><span class="p">,</span> <span class="n">X</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="n">X</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="n">X</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="n">s0</span><span class="p">]</span>
    
    <span class="c1"># Add in the state space for the dispersed phase particles</span>
    <span class="n">q</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">dispersed_phases</span><span class="o">.</span><span class="n">particles_state_space</span><span class="p">(</span><span class="n">particles</span><span class="p">,</span> <span class="n">nbe</span><span class="p">))</span>
    
    <span class="c1"># Add the ambient concentrations of the dispersed-phase chemicals</span>
    <span class="n">ca</span> <span class="o">=</span> <span class="n">profile</span><span class="o">.</span><span class="n">get_values</span><span class="p">(</span><span class="n">X</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="n">chem_names</span><span class="p">)</span>
    <span class="n">q</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">Mj</span> <span class="o">/</span> <span class="n">rho_j</span> <span class="o">*</span> <span class="n">ca</span><span class="p">)</span>
    
    <span class="c1"># Add in the tracers discharged with the jet</span>
    <span class="n">q</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">Mj</span><span class="o">*</span><span class="n">cj</span><span class="p">)</span>
    
    <span class="c1"># Return the complete initial conditions</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">t</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">(</span><span class="n">q</span><span class="p">))</span></div>


<div class="viewcode-block" id="zfe_volume_flux"><a class="viewcode-back" href="../autodoc/bpm/lmp.zfe_volume_flux.html#lmp.zfe_volume_flux">[docs]</a><span class="k">def</span> <span class="nf">zfe_volume_flux</span><span class="p">(</span><span class="n">profile</span><span class="p">,</span> <span class="n">X0</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">Vj</span><span class="p">,</span> <span class="n">Sj</span><span class="p">,</span> <span class="n">Tj</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute the volume flux of continous phase discharge fluid at the release</span>
<span class="sd">    </span>
<span class="sd">    If the release includes continous phase fluid, this function computes</span>
<span class="sd">    the flow rate and geometry of the release</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    profile : `ambient.Profile` object</span>
<span class="sd">        The ambient CTD object used by the single bubble model simulation.</span>
<span class="sd">    X0 : ndarray</span>
<span class="sd">        Release location (x, y, z) in (m)</span>
<span class="sd">    R : float</span>
<span class="sd">        Radius for the equivalent circular cross-section of the release </span>
<span class="sd">        (m)</span>
<span class="sd">    Vj : float</span>
<span class="sd">        Scalar value of the magnitude of the discharge velocity for </span>
<span class="sd">        continuous phase fluid in the discharge.  This variable should be </span>
<span class="sd">        0 or None for a pure multiphase discharge.</span>
<span class="sd">    Sj : float</span>
<span class="sd">        Salinity of the continuous phase fluid in the discharge (psu)</span>
<span class="sd">    Tj : float</span>
<span class="sd">        Temperature of the continuous phase fluid in the discharge (T)</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    Q : Volume flux of continuous phase fluid at the discharge (m^3/s)</span>
<span class="sd">    A : Cross-sectional area of the discharge (M^2)</span>
<span class="sd">    X : ndarray</span>
<span class="sd">        Release location (x, y, z) in (m)</span>
<span class="sd">    Tj : float</span>
<span class="sd">        Temperature of the continuous phase fluid in the discharge (T)</span>
<span class="sd">    Sj : float</span>
<span class="sd">        Salinity of the continuous phase fluid in the discharge (psu)</span>
<span class="sd">    Pj : float</span>
<span class="sd">        Pressure at the discharge (Pa)</span>
<span class="sd">    rho_j : float</span>
<span class="sd">        Density of the continous phase fluid in the discharge (kg/m^3)</span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># The Lagrangian plume model starts at the discharge.</span>
    <span class="n">X</span> <span class="o">=</span> <span class="n">X0</span>
    
    <span class="c1"># Get the jet density from the discharge characteristics</span>
    <span class="n">Ta</span><span class="p">,</span> <span class="n">Sa</span><span class="p">,</span> <span class="n">P</span> <span class="o">=</span> <span class="n">profile</span><span class="o">.</span><span class="n">get_values</span><span class="p">(</span><span class="n">X</span><span class="p">[</span><span class="mi">2</span><span class="p">],</span> <span class="p">[</span><span class="s1">&#39;temperature&#39;</span><span class="p">,</span> <span class="s1">&#39;salinity&#39;</span><span class="p">,</span> 
        <span class="s1">&#39;pressure&#39;</span><span class="p">])</span>
    <span class="n">rho_j</span> <span class="o">=</span> <span class="n">seawater</span><span class="o">.</span><span class="n">density</span><span class="p">(</span><span class="n">Ta</span><span class="p">,</span> <span class="n">Sa</span><span class="p">,</span> <span class="n">P</span><span class="p">)</span>
    
    <span class="c1"># Pressure at the discharge is the ambient pressure</span>
    <span class="n">Pj</span> <span class="o">=</span> <span class="n">P</span>
    
    <span class="c1"># The discharge area if the full port area</span>
    <span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">pi</span> <span class="o">*</span> <span class="n">R</span><span class="o">**</span><span class="mi">2</span>
    
    <span class="c1"># Compute the volume flux of discharge fluid</span>
    <span class="n">Q</span> <span class="o">=</span> <span class="n">A</span> <span class="o">*</span> <span class="n">Vj</span>
    
    <span class="c1"># Return the initial conditions with salinity and temperature of the </span>
    <span class="c1"># discharge equal to the jet values</span>
    <span class="k">return</span> <span class="p">(</span><span class="n">Q</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">X</span><span class="p">,</span> <span class="n">Tj</span><span class="p">,</span> <span class="n">Sj</span><span class="p">,</span> <span class="n">Pj</span><span class="p">,</span> <span class="n">rho_j</span><span class="p">)</span></div>


<div class="viewcode-block" id="flux_to_velocity"><a class="viewcode-back" href="../autodoc/bpm/lmp.flux_to_velocity.html#lmp.flux_to_velocity">[docs]</a><span class="k">def</span> <span class="nf">flux_to_velocity</span><span class="p">(</span><span class="n">Q</span><span class="p">,</span> <span class="n">A</span><span class="p">,</span> <span class="n">phi</span><span class="p">,</span> <span class="n">theta</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Convert fluid flow rate to three-component velocity </span>
<span class="sd">    </span>
<span class="sd">    Computes the three-component velocity (u, v, w) along the Cartesian </span>
<span class="sd">    directions (x, y, depth) from the flow rate, cross-sectional area, and the</span>
<span class="sd">    orientation (phi and theta).</span>
<span class="sd">    </span>
<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    Q : Volume flux of continuous phase fluid at the discharge (m^3/s)</span>
<span class="sd">    A : Cross-sectional area of the discharge (M^2)</span>
<span class="sd">    phi : float</span>
<span class="sd">        Vertical angle from the horizontal for the discharge orientation </span>
<span class="sd">        (rad in range +/- pi/2)</span>
<span class="sd">    theta : float</span>
<span class="sd">        Horizontal angle from the x-axis for the discharge orientation.  </span>
<span class="sd">        The x-axis is taken in the direction of the ambient current.  </span>
<span class="sd">        (rad in range 0 to 2 pi)</span>
<span class="sd">    </span>
<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    Uj : ndarray</span>
<span class="sd">        Vector of the three-component velocity of continous phase fluid in </span>
<span class="sd">        the jet (u, v, w) in the Cartesian direction (x, y, depth) </span>
<span class="sd">    </span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="c1"># Get the velocity along the jet centerline</span>
    <span class="n">Vj</span> <span class="o">=</span> <span class="n">Q</span> <span class="o">/</span> <span class="n">A</span>
    
    <span class="c1"># Project jet velocity on the three component directions (i, j, k)</span>
    <span class="n">Uj</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span>
    <span class="n">Uj</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">phi</span><span class="p">)</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">theta</span><span class="p">)</span> <span class="o">*</span> <span class="n">Vj</span>
    <span class="n">Uj</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">phi</span><span class="p">)</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">theta</span><span class="p">)</span> <span class="o">*</span> <span class="n">Vj</span>
    <span class="n">Uj</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">phi</span><span class="p">)</span> <span class="o">*</span> <span class="n">Vj</span>
    
    <span class="c1"># Return the velocity vector</span>
    <span class="k">return</span> <span class="n">Uj</span></div>

</pre></div>

          </div>
          
        </div>
      </div>
      <div class="sphinxsidebar" role="navigation" aria-label="main navigation">
        <div class="sphinxsidebarwrapper">
<h1 class="logo"><a href="../index.html">Texas A&M Oil spill / Outfall Calculator</a></h1>








<h3>Navigation</h3>
<ul>
<li class="toctree-l1"><a class="reference internal" href="../user_manual.html">TAMOC User Manual</a></li>
<li class="toctree-l1"><a class="reference internal" href="../unit_tests.html">Unit Tests</a></li>
<li class="toctree-l1"><a class="reference internal" href="../glossary.html">Glossary</a></li>
</ul>
<ul>
<li class="toctree-l1"><a class="reference internal" href="../bugs.html">Reporting Bugs</a></li>
<li class="toctree-l1"><a class="reference internal" href="../readme.html">Read Me File</a></li>
<li class="toctree-l1"><a class="reference internal" href="../release.html">Release Notes</a></li>
<li class="toctree-l1"><a class="reference internal" href="../license.html">License</a></li>
</ul>

<div class="relations">
<h3>Related Topics</h3>
<ul>
  <li><a href="../index.html">Documentation overview</a><ul>
  <li><a href="index.html">Module code</a><ul>
  </ul></li>
  </ul></li>
</ul>
</div>
<div id="searchbox" style="display: none" role="search">
  <h3 id="searchlabel">Quick search</h3>
    <div class="searchformwrapper">
    <form class="search" action="../search.html" method="get">
      <input type="text" name="q" aria-labelledby="searchlabel" />
      <input type="submit" value="Go" />
    </form>
    </div>
</div>
<script>$('#searchbox').show(0);</script>








        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="footer">
      &copy;2020, Scott A. Socolofsky.
      
      |
      Powered by <a href="http://sphinx-doc.org/">Sphinx 2.4.4</a>
      &amp; <a href="https://github.com/bitprophet/alabaster">Alabaster 0.7.12</a>
      
    </div>

    

    
  </body>
</html>